An introduction to open-source R software for wearable device data processing and analysis

Introduction

R is a functional programming language for statistical computing and graphics.¹
RStudio is an integrated development environment for R
There are >20,000 available packages on CRAN extending the utility of R
This presentation was created in Quarto using R

R is a functional programming language for statistical computing and graphics. A functional programming language means that you can extend the functionality of the programming language by writing functions designed to accomplish specific tasks. RStudio is an integrated development environment for R. However, it is not the only way to interact with R on your computer. For example, some users prefer to interact with R in a separate R console or integrated development environment like VS Code. R can also be run from your terminal. Common functions are often assembled into a package for repeated use and there are > 20,000 packages on CRAN that extend the utility of R with even more packages available on GitHub. This presentation will introduce you to open-source R software for wearable device data processing and analysis. The presentation was assembled in Quarto which is an open source technical publishing system to create articles, websites, blogs, books, and presentations.

Why R?

Many open-source physical activity packages available on CRAN or GitHub
Open-source code free to view on GitHub
Reticulate is a Python-R interface that makes Python packages available to R users
Customized scripts or packages can meet the needs of a research study.

There are many open-source physical activity packages available on CRAN or GitHub. Some commonly used packages include read.gt3x for reading binary data, PhysicalActivity for wear time processing, and GGIR for analyzing raw acceleration sleep and physical activity data. The compositions package offers an easy way to analyze movement data as a composition and I’ve been working on developing agcounts and iFitbit. The agcounts package allows users to calculate ActiGraph counts from raw acceleration data and iFitbit interacts with the Fitbit API to extract data from research participants. Open-source code from these packages are free to view on GitHub and I’ve found that it is a great way to learn. Issues on GitHub are like discussion forums where users can discuss problems with the code and request new features. A GitHub pull request allows users to suggest, discuss, and review changes to the code and version control provides a method for tracking the changes over time. Additionally, R can be extended to Python packages for wearable device data processing and analysis through a Python-R interface called reticulate. R also allows for easy customization of scripts or packages to meet the needs of a research study.

Objectives

The purpose of today’s presentation is to:

Introduce you to R packages used in physical activity research
Provide examples with code for how physical activity data can be processed, analyzed, and visualized in R.

Getting Started

Downloading R and RStudio https://posit.co/download/rstudio-desktop/
The devtools package allows you to install packages from GitHub

# Install the devtools package from CRAN
install.packages("devtools")

Download the Presentation

Install this presentation into RStudio as an R package

# Install the presentationsR package
options(timeout = 999999) # Increase to avoid timeout issues
devtools::install_github("bhelsel/presentationsR")
presentationsR::render_presentation(outdir = "/full/path/to/folder", name = "csacsm24")

Visit GitHub to download as a compressed .zip file

https://github.com/bhelsel/presentationsR

Wearable Devices

Accelerometry was first introduced in the 1980’s²
Over the past ~40 years:
- The physical size of portable accelerometers is smaller
- More onboard data storage capacity
- Refined data collection and processing protocols
- Open-source software is becoming increasingly prevalent

Open-source Software

Has an important role in advancing scientific reproducibility and rigor
Can be developed in a collaborative public fashion
Software needs to be findable, accessible, interoperable, and reusable
Many scientists who develop or use open-source software are self-taught
Comprehensive, didactic documentation and training opportunities are important

A Basic Accelerometry Workflow

Read in the acceleration data from ActiGraph .gt3x and .agd files
Convert the raw acceleration data to ActiGraph counts
Identify non-wear time using the Choi algorithm
Apply cut-points and separate acceleration data into intensity bands³
Create compositional variables for analysis and visualizations

We are going to look at a basic accelerometry workflow in R. This workflow may look different for you depending on the type of monitor you use or the metrics that you are interested in from the accelerometer. First, I’ll show you how to read in the acceleration data from the ActiGraph .gt3x and .agd files. The .gt3x files contain binary data for compact storage and the .agd files are SQL databases. Second, I’ll show you how to convert raw acceleration data to ActiGraph counts using ActiGraph’s open-source count algorithm. The count algorithm was proprietary up until 2022 when ActiGraph made it publicly available in Python and our team translated it into an R package. Third, we will apply the Choi algorithm to detect non-wear time, apply cut-points, and separate the acceleration data into intensity bands. Intensity bands are an alternative to cut-points that separates the acceleration data into predetermined bins. Finally, we will talk about analyzing and visualizing accelerometry data as a compositional variable. Originally, this was an approach for analyzing mineral compositions in chemistry, but is now being used to study time-use behaviors since sleep, sedentary time, and physical activity make up an entire day. This is a powerful analytic approach that allows you study the effect of time-use reallocation on health outcomes.

Read in the raw acceleration data

# Binary .gt3x files can be read with the read.gt3x pacakge
gt3xFile <- system.file("extdata/example3min.gt3x", package = "presentationsR")
gt3xData <- read.gt3x::read.gt3x(gt3xFile, asDataFrame = TRUE)
head(gt3xData)

Sampling Rate: 100Hz
Firmware Version: 1.7.2
Serial Number Prefix: TAS
                 time      X     Y     Z
1 2023-06-13 08:34:00 -0.020 0.008 1.031
2 2023-06-13 08:34:00  0.000 0.004 1.035
3 2023-06-13 08:34:00  0.008 0.004 1.031
4 2023-06-13 08:34:00  0.008 0.004 1.035
5 2023-06-13 08:34:00  0.008 0.004 1.031
6 2023-06-13 08:34:00  0.016 0.008 1.027

sprintf("Data has %s rows and %s columns", nrow(gt3xData), ncol(gt3xData))

[1] "Data has 18000 rows and 4 columns"

One of the easiest and quickest way to read in the binary data to R is through the read.gt3x R package. This package integrates C++ code into the reading functions which provides increased data processing speed. presentationsR has example .gt3x and .agd data contained within the R package which can be obtained by the system.file function when installed. The system.file function retrieves the path name for the file’s location which can also be copied manually if you download the compressed .zip file containing the contents of the R package. Next, you can read in the raw acceleration data to your R environment using the read.gt3x function from the read.gt3x package. The two colons attach the function to the package. Alternatively, R packages can be loaded with the library function. The asDataFrame argument converts the data set into a data.frame object rather than a matrix which simplifies data processing. The data has 18,000 rows and contains a time variable as well as acceleration data from the X, Y, and Z axes.

Example Output

Using Reticulate

ActiGraph has the pygt3x file reader in a Python package https://github.com/actigraph/pygt3x which has subtle differences from the read.gt3x R package.

if(reticulate::py_module_available("pygt3x")){
  `%as%` <- reticulate::`%as%`
  Reader <- reticulate::import("pygt3x.reader", convert = FALSE)
  pd <- reticulate::import("pandas")
  reset_index <- pd$core$frame$DataFrame$reset_index
  FileReader <- Reader$FileReader
  to_pandas <- FileReader$to_pandas
  with(FileReader(gt3xFile) %as% reader, {
    gt3xDataPy = reset_index(to_pandas(reader))
  })
  colnames(gt3xDataPy) <- c("time", "X", "Y", "Z")
  gt3xDataPy$time <- as.POSIXct(
    gt3xDataPy$time, origin = "1970-01-01 00:00:00", tz = "UTC"
    )
}

An additional option for reading in data from the .gt3x files is by using ActiGraph’s pygt3x reader. This reader offers a flexible solution to reading data from the .gt3x files and aligns with ActiGraph’s recommendations. It can also calibrate the data from the CentrePoint Insight Watch and handle the .agdc file types. However, this reader is written as a Python package which requires users to work in both programming languages. The reticulate package offers a Python-R interface that can be used to call Python packages in R. It does require users to navigate Python installation and environment setup before using reticulate to access Python packages, but here is an example of what that might look like in R after the set up is completed. First, you’d make sure that the pygt3x module is available. The import function loads the necessary Python packages and classes before attaching them to an R object with the $ sign operator as seen in lines 5-7. The data can then be read by FileReader as displayed in lines 8-10. Accessing Python packages via reticulate greatly expands the utility of R for processing wearable device data.

Converting to ActiGraph Counts

ActiGraph’s count algorithm was made publicly available in 2022.⁴
The agcounts Python module was converted to an R package in 2022 and is available on GitHub https://github.com/bhelsel/agcounts.

# ::: allows you to call hidden and non-exported functions
sf <- agcounts:::.get_frequency(gt3xData)
resampled_data <- agcounts:::.resample(gt3xData, frequency = sf)
filtered_data <- agcounts:::.bpf_filter(resampled_data)
trimmed_data <- agcounts:::.trim_data(filtered_data, lfe_select = FALSE)
resampled_10Hz_data <- agcounts:::.resample_10hz(trimmed_data)
epoch_counts <- agcounts:::.sum_counts(resampled_data, epoch = 60)

Once we’ve used a reader to import the raw acceleration data into R, we can convert the data to ActiGraph counts using the agcounts package. The steps of the ActiGraph count algorithm are shown on this slide. To start, we need to know the sampling frequency. If you don’t have this recorded as part of your study information, we’ve created a method for retrieving the sampling frequency from within the agcounts R package using the .get_frequency function. It is not an exported function from the package since it’s mostly for internal use, but using 3 colons will allow you to call hidden and non-exported functions. If you pass in the raw acceleration data to the .get_frequency function, you will see that the sampling frequency is 100 Hz. Then, you can execute each step of ActiGraph’s count algorithm to down sample the data to 30 Hz, pass the data through a band pass filter, add a normal or low frequency extension filter, resample the filtered data to 10 Hz, and add the counts over a user-specified epoch such as 60 seconds.

The Power of Functions

.bpf_filter <- function(downsample_data, verbose = FALSE){
  if(verbose){
    print("Filtering Data")
  }
  a = as.numeric(.coefficients$output_coefficients)
  b = as.numeric(.coefficients$input_coefficients)
  zi <- gsignal::filter_zi(filt = b, a = a)
  zi <- matrix(rep(zi, 3), nrow = 3, byrow = 3) * downsample_data[, 1]
  rownames(zi) <- c("X", "Y", "Z")
  bpf_data_x <- gsignal::filter(filt = b, a = a, x = downsample_data["X", ], zi = zi["X", ])
  bpf_data_y <- gsignal::filter(filt = b, a = a, x = downsample_data["Y", ], zi = zi["Y", ])
  bpf_data_z <- gsignal::filter(filt = b, a = a, x = downsample_data["Z", ], zi = zi["Z", ])
  bpf_data <- matrix(rbind(bpf_data_x$y, bpf_data_y$y, bpf_data_z$y), nrow = 3)
  bpf_data = ((3.0 / 4096.0) / (2.6 / 256.0) * 237.5) * bpf_data # 17.127404 is used in ActiLife and 17.128125 is used in firmware.
  rownames(bpf_data) <- c("X", "Y", "Z")
  return(bpf_data)
}

A Wrapper Function

A function with a main purpose of calling secondary functions or code to ease implementation.

agcounts::get_counts(gt3xFile, epoch = 30, lfe_select = TRUE, parser = "read.gt3x")

                 time Axis1 Axis2 Axis3 Vector.Magnitude
1 2023-06-13 08:34:00  1723  1384  2206             3123
2 2023-06-13 08:34:30  1032  1842  1455             2564
3 2023-06-13 08:35:00  1055  2074  1420             2726
4 2023-06-13 08:35:30   871  1974  1520             2639
5 2023-06-13 08:36:00   924  1717  1278             2331
6 2023-06-13 08:36:30  1377  1764  1480             2683

agcounts::get_counts(gt3xFile, epoch = 60, lfe_select = FALSE, parser = "read.gt3x")

                 time Axis1 Axis2 Axis3 Vector.Magnitude
1 2023-06-13 08:34:00  2606  3116  3542             5389
2 2023-06-13 08:35:00  1738  3943  2840             5161
3 2023-06-13 08:36:00  2172  3363  2646             4799

SQL Database

ActiGraph’s .agd files are SQL databases
Data from a SQL database is loaded with the DBI package

agdFile <- system.file("extdata/example20hours1sec.agd", package = "presentationsR")
# Open the SQL database connection
con <- DBI::dbConnect(RSQLite::SQLite(), agdFile)
# Read the data into R from the SQL database
agdData <- DBI::dbReadTable(con, "data")
# Close the SQL database
DBI::dbDisconnect(con)

                 time axis1 axis2 axis3   vm
1 2023-04-24 15:02:00  1235  1646  1550 3079
2 2023-04-24 15:03:00   808   931  1102 1846
3 2023-04-24 15:04:00   361   410  1437 1628
4 2023-04-24 15:05:00  1060  1004  1927 2639
5 2023-04-24 15:06:00   697  1079   755 1700
6 2023-04-24 15:07:00   126   249   482  585

Wear Time Marking

Choi et al. (2011)⁵

agdData <- 
  agdData |>
  PhysicalActivity::wearingMarking(
    frame = 90,
    perMinuteCts = 1, 
    TS = "time",
    cts = "vm",
    allowanceFrame = 2,
    newcolname = "wear"
  )
head(agdData)

                 time axis1 axis2 axis3   vm wear weekday days
1 2023-04-24 15:02:00  1235  1646  1550 3079    w  Monday    1
2 2023-04-24 15:03:00   808   931  1102 1846    w  Monday    1
3 2023-04-24 15:04:00   361   410  1437 1628    w  Monday    1
4 2023-04-24 15:05:00  1060  1004  1927 2639    w  Monday    1
5 2023-04-24 15:06:00   697  1079   755 1700    w  Monday    1
6 2023-04-24 15:07:00   126   249   482  585    w  Monday    1

Wear time marking can be completed using the PhysicalActivity R package. It was written by Dr. Leena Choi and allows users to classify wear and non-wear time for accelerometer data on an epoch-by-epoch basis. There are several arguments within the function that can be adjusted based on the user’s preferences. These include the frame, perMinuteCts, and allowanceFrame. The frame is the time interval to be considered and allowanceFrame is the size of time interval that zero counts are allowed. These default to the Choi et al. (2011) defaults of 90 minutes and 2 minutes, respectively. perMinuteCts are the number of data rows per minute which defaults to 1-second epochs or 60 rows per minute. Here we have a value of 1 because our data is currently using 1-minute epochs. The other arguments within the wearingMarking function are related to column names. TS identifies the column name of the time stamped variable, cts is the column containing the counts that you want to use to calculate non-wear, and newcolname is the name of the column in the new data set that will contain the wear and non-wear information.

Cut-points and Intensity Bands

Montoye et al (2020)⁶

classified_data <- 
  agdData |>
  dplyr::mutate(intensity = cut(
    vm, breaks = c(0, 2860, 3941, Inf), # Montoye Wrist
    labels = c("Sedentary", "Light", "MVPA"),
    include.lowest = TRUE, right = FALSE)) |>
  dplyr::mutate(intensity_bands = cut(
    vm, breaks = seq(0, 10500, 1500),
    labels = c("0_1500", "1501_3000", "3001_4500", "4501_6000", 
               "6001_7500", "7501_9000", "> 9000"),
    include.lowest = TRUE, right = FALSE))

head(classified_data[, -c(2:4)])

                 time   vm wear weekday days intensity intensity_bands
1 2023-04-24 15:02:00 3079    w  Monday    1     Light       3001_4500
2 2023-04-24 15:03:00 1846    w  Monday    1 Sedentary       1501_3000
3 2023-04-24 15:04:00 1628    w  Monday    1 Sedentary       1501_3000
4 2023-04-24 15:05:00 2639    w  Monday    1 Sedentary       1501_3000
5 2023-04-24 15:06:00 1700    w  Monday    1 Sedentary       1501_3000
6 2023-04-24 15:07:00  585    w  Monday    1 Sedentary          0_1500

Aggregating Intensity Data

intensity_binary <- cbind(
  classified_data, 
  sjmisc::to_dummy(classified_data$intensity)
  )

agdDataHour <- 
  intensity_binary |>
  dplyr::group_by(time = format(time, "%Y-%m-%d %H")) |>
  dplyr::summarise_at(dplyr::vars(intensity_1:intensity_3), sum) |>
  `colnames<-`(c("time", "Sedentary", "Light", "MVPA"))

head(agdDataHour)

# A tibble: 6 × 4
  time          Sedentary Light  MVPA
  <chr>             <dbl> <dbl> <dbl>
1 2023-04-24 15        48     5     5
2 2023-04-24 16        51     5     4
3 2023-04-24 17        27     7    26
4 2023-04-24 18        34    17     9
5 2023-04-24 19        49     5     6
6 2023-04-24 20        35     9    16

Now that we have our classified data, it is time to aggregate the intensity data so we can start looking at summaries. Usually, we’d consider aggregating the data to a day-level estimate, but the sample .agd file that we loaded only has ~20 hours of data. Thus, as an example, we will aggregate the data to an hour-level estimate of movement. I’m showing how to do this with our cut-point data but a similar approach could be used for the intensity bands. The to_dummy function in the sjmisc R package can take a character or factor variable and convert it to a set of dummy variables based on the variable’s unique values. These columns can be added to the original data set using cbind. Next, the user can use the group_by and summarise_at functions from the dplyr package to group and sum the data by hour.

Compositional Data

# Change to composition data
comp <- compositions::acomp(agdDataHour[, 2:4])
head(round(comp, 3))

     Sedentary Light  MVPA
[1,]     0.828 0.086 0.086
[2,]     0.850 0.083 0.067
[3,]     0.450 0.117 0.433
[4,]     0.567 0.283 0.150
[5,]     0.817 0.083 0.100
[6,]     0.583 0.150 0.267

# Calculate geometric means using 60 minutes as the total
unname(compositions::clo(mean(comp, robust = FALSE), total = 60))

[1] 50.667662  4.170343  5.161995

Analysis: Isometric Log Ratio (ILR) Transformations

Van den Boogaart et al. (2008) Compositions: A unified R package to analyze compositional data⁷

# Sequential Binary Partition
sbp <- matrix(
  c(-1, -1,
    -1, +1,
    +1,  0),
  byrow = TRUE,
  ncol = 2)
# Base of the centered log ratio (CLR) plane
psi <- compositions::gsi.buildilrBase(sbp)
# Isometric log-ratio (ILR) transformed compositions
ilr_transformed <- compositions::ilr(comp, V = psi)
head(ilr_transformed, 5)

           [,1]       [,2]
[1,] -0.9233609 -1.5993080
[2,] -1.1303068 -1.6421761
[3,]  0.5202904 -0.9545423
[4,] -0.8022588 -0.4901291
[5,] -0.7829138 -1.6138881

We can analyze compositional data in regressions, but including all of the compositional variables introduces multicollinearity into the model since the time-use behaviors will be highly correlated. We can use an isometric log-ratio transformation to map a D-part composition onto a D-1 dimensional euclidean vector to be analyzed by all of the classical multivariate tools. To accomplish this using the compositions package, we need to first write a sequential binary partition containing +1, 0, and -1 values according to our analysis plans. A positive 1 indicates that the variable should be included in the numerator, a value of 0 excludes the variable, and a negative 1 adds the variable to the denominator. The sequential binary partition is arbitrary and can be adjusted based on your research question. For example, if I decided that I want to look at MVPA relative to sedentary time, I could change the -1 to 0 in the second row and first column. Similarly, I could change row 3, column 2 to +1 to look at physical activity (light and MVPA) relative to sedentary time. The math is straightforward but involves some understanding of log transformations and matrix multiplication. I added a few basic functions in this package to simplify the code used in the compositions R package to help improve your understanding of the isometric log ratio transformation.

Behind the Scenes: The ILR Function

isPos = (sbp > 0)
isNeg = (sbp < 0)
m <- matrix(1, nrow(sbp), nrow(sbp))
nPos = m %*% isPos
nNeg = m %*% isNeg
sbp_t = (isPos * nNeg - isNeg * nPos)
nn = sapply(1:ncol(sbp_t), function(i) {
  1/sqrt(sbp_t[, i] %*% sbp_t[, i])
})
nn = matrix(nn, ncol = ncol(sbp_t), nrow = nrow(sbp_t), byrow = TRUE)
psi = sbp_t * nn
LOG <- unclass(log(ifelse(comp > 0, comp, 1)))
clr <- ifelse(comp > 0, LOG - rowSums(LOG)/rowSums(comp > 0), 0)
ilr <- clr %*% psi
all(ilr == ilr_transformed)

[1] TRUE

Analysis: An Alternative ILR Transformation

\[ z1 = \sqrt{\frac{2}{3}} \cdot \ln\left(\frac{MVPA}{\sqrt[2]{SB \cdot LPA}}\right) \]

\[ z2 = \sqrt{\frac{1}{2}} \cdot \ln\left(\frac{LPA}{\sqrt[1]{SB}}\right) \]

Templ et al. (2011) robCompositions: An R-package for Robust Statistical Analysis of Compositional Data⁸

X <- comp[, c("MVPA", "Light", "Sedentary")]
utils::head(robCompositions::pivotCoord(X))

  MVPA_Li-Se   Light_Se
1 -0.9969690 -1.5993080
2 -1.1866392 -1.6421761
3  0.0565334 -0.9545423
4 -0.9349550 -0.4901291
5 -0.8689403 -1.6138881
6 -0.3379092 -0.9603383

all.equal(MoveKC::ilr_pivotCoord(X), robCompositions::pivotCoord(X))

[1] TRUE

Next Steps

Visualize the data using ternary plots
Use different multivariate tools
Raw acceleration data
Explore isotemporal substitution models to determine the effects of time-use reallocation

Now that you’ve completed the data processing, you can start to explore the composition using visualizations and multivariate tools. A plot that is useful to visualize compositional data is a ternary plot similar to the one pictured on this slide. A ternary plot is a graphical representation of data consisting of 3 components. In this example, our components include sleep, sedentary time, and physical activity with an average composition of 32.2% of the day in sleep, 56.6% of the day spent doing sedentary activities, and 11.2% of the day spent doing light or moderate-to-vigorous physical activity. You can also start to explore raw acceleration data and isotemporal substitution models which look at how the reallocation of time between sleep, physical activity, and sedentary time might affect your health outcomes. This may help answer causal questions in observational studies about dose-effect responses to increasing physical activity.

Fitbit API

Fitbit provides access to user data via their application programming interface (API)
Third-party companies exist to make access to the Fitbit API easier for investigators
- Pilot studies lack the budget to pay for the service
- Participant data is linked to another platform
Tools in R make it possible to retrieve Fitbit data and it can be automated!

iFibit

We created iFitbit to simplify working with the Fitbit API
It allows you to extract activity, exercise, sleep, and heart rate data to:
- Store in a SQL database
- Run interactive reports

To install:

devtools::install_github("bhelsel/iFitbit", build_vignettes = TRUE)
vignette("fitbit-application", package = "iFitbit")

We created iFitbit to simplify working with the Fitbit API. It allows you to extract activity, exercise, sleep, and heart rate data to store in a SQL database and run interactive reports for the research team. You can install it using the devtools::install_github function that we used earlier in this presentation. However, if you want to install the “getting started” vignette, you will want to add build_vignettes = TRUE into the function. This will provide access to a vignette that will walk you through some of the steps for registering as a Fitbit developer. The README.md page on GitHub has information on interacting with the R package and code examples, so I will not do a deep dive into the package here. Instead, I will focus on a workflow for extracting, processing, and automating the Fitbit reports.

Example Data from the Fitbit API

https://dev.fitbit.com/build/reference/web-api/

fb_file <- system.file("extdata/exampleFitbit.db", package = "presentationsR")
con <- DBI::dbConnect(RSQLite::SQLite(), fb_file)
tbl_n <- DBI::dbListTables(con) # List table names
data <- lapply(tbl_n, function(x) DBI::dbReadTable(con, name = x))
names(data) <- tbl_n
tbl_n

[1] "activities"   "device"       "exercise_log" "heart"        "sleep"

Included in the CASCSM24 presentation is an example database containing data extracted from the Fitbit Web API. It has activity data, information about the device, an exercise log, intensity by heart rate, and sleep data. This is the information we have used in previous studies, but is a lot of other data that can be extracted from the API like body weight, nutrition, heart rate variability, fitness, social network, and more. A complete list of data available can be found by navigating to the link posted on this slide. An example database can be loaded into your R environment by using the DBI::dbConnect function that we used earlier for the ActiGraph .agd files. The lapply function here loops through the tables in the SQL database and returns them as a list to your R environment. This is an easy way to load all the data at the same time.

Activities

date	distance	elevation	floors	minutesSedentary	minutesLightlyActive	minutesFairlyActive	minutesVeryActive	steps
2022-02-28	1.840	6	6	1169	206	45	20	9412
2022-03-01	2.397	22	4	1319	86	26	9	8150
2022-03-02	2.028	12	6	1240	161	34	5	8496
2022-03-03	2.565	10	1	1249	124	52	15	6076
2022-03-04	1.197	28	4	1192	206	25	17	10728
2022-03-05	0.820	10	3	1190	212	21	17	9668
2022-03-06	2.392	5	2	1226	175	35	4	8579
2022-03-07	2.382	4	0	1256	161	18	5	8543
2022-03-08	0.733	4	2	1257	147	24	12	7326
2022-03-09	1.245	33	6	1297	113	20	10	5612

Exercise Log

date	time	type	duration	steps	calories	hr
2022-02-28	20:45:15	Swim	29.22	8	527	80
2022-03-01	02:35:04	Swim	12.41	15	767	123
2022-03-02	02:19:55	Swim	35.55	8	531	91
2022-03-03	00:20:32	Sport	15.21	334	669	99
2022-03-04	16:11:54	Walk	11.26	734	494	94
2022-03-05	03:21:53	Yoga	28.87	28	413	121
2022-03-06	22:04:07	Strength	13.84	321	554	112
2022-03-07	03:49:23	Workout	10.49	421	568	98
2022-03-08	10:26:17	Aerobic Workout	26.82	644	555	86
2022-03-09	19:34:00	Sport	29.74	148	608	115

Sleep and Intensity by Heart Rate

date	sleep	nonwear	sedentary	very.light	light	moderate	vigorous	near.maximal	mvpa
2022-02-28	397	563	519	174	152	22	4	6	32
2022-03-01	401	620	507	138	137	22	11	5	38
2022-03-02	447	613	518	153	120	22	8	6	36
2022-03-03	391	648	457	181	120	16	15	3	34
2022-03-04	421	653	410	209	124	25	14	5	44
2022-03-05	359	485	588	201	128	23	10	5	38
2022-03-06	445	576	469	184	149	50	7	5	62
2022-03-07	407	571	541	152	140	26	8	2	36
2022-03-08	399	543	586	157	111	37	4	2	43
2022-03-09	372	557	491	205	131	36	14	6	56

Heart Rate Intraday Function

get_fitbit_heart_intraday(
  token.pathname,
  start.date = "2022-02-28",
  end.date = "2022-03-20",
  detail.level = "1min",
  age = 25,
  max_hr_formula = "220 - age",
  intensity = c(0, 0.2, 0.3, 0.4, 0.6, 0.9, 1),
  intensity_labels = c("sedentary", "very.light", "light", "moderate", "vigorous", "near.maximal"),
  heart_rate_method = 2,
  rest.hr = 60,
  do.parallel = TRUE,
  returnData = FALSE,
  returnRawData = FALSE,
  toSQL = TRUE,
  verbose = TRUE
)

Here is the heart rate intraday function within the iFitbit package to give you an idea of some of the customization that can happen. The token.pathname is the access token to the user’s account that will be saved to your computer as part of the authentication process. The start.date and end.date are the range of dates for which you want the heart rate data extracted and detail.level is the epoch level desired which can be 1 second, 1 minute, 5 minutes, or 15 minutes. Lines 6-11 are where the customization of the heart rate intensity classification will happen. Users can change the age-based heart rate formula or alter any of the intensities or their labels. The heart_rate_method allows users to easily switch between heart rate reserve or maximal heart rate formulas and resting heart rate can be inputed into the function or extracted from the Fitbit data if it is available. The do.parallel argument helps speed up processing time by using n-1 cores on your computer to process the intraday data since there could be up to 86,400 observations per day if the heart rate data are extracted in 1-second level epochs. Finally, the user can choose to return the raw or summarized data to their R environment or write the data to a SQL database assuming one was created during the authentication process.

Intraday Visualization of Heart Rate Data

Interday Visualization of Heart Rate Data

Example Report

Automating the Process

cmd <- cronR::cron_rscript("pathname/to/script.R", log_append = FALSE)
cronR::cron_add(
  command = cmd,
  frequency = '0 9 * * 1-5', # 9 am Monday-Friday
  id = 'iFitbitReport',
  user = "username",
  description = 'Extracts data from the Fitbit API and runs interactive reports',
  tags = c('fitbit', 'interactive', 'exercise', 'sleep'))

cronR::cron_ls("iFitbitReport") # List cron jobs
cronR::cron_rm('iFitbitReport') # Remove cron job

Conclusions

Many scientists view software as an important part of their research
Scientists developing software may benefit from open-source resources like GitHub to:
- Enhances the learning process by allowing others to inspect, modify, and enhance their code
- Increases reproducibility and drives innovation
Customized software for wearable devices can save time and increase collaboration in physical behavior research

In conclusion, many scientists view software as an important part of their research. I hope the presentation today demonstrated some of the value that open-source software can have in research using wearable devices. Scientists who are developing software, even at a beginner level, may benefit from open-source resources like GitHub that allow users to track and discuss changes to the software developed. This may enhance the learning process for the scientist developing the software by allowing others inspect, modify, and enhance their code. This also increases reproducibility and drives innovation. Finally, customized software for wearable devices can save time and increase collaboration in physical behavior research. It has been a skillset that I developed as a postdoctoral fellow at the University of Kansas and it has led to collaborations with investigators at various institutions.

References

1.

R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing; 2021. https://www.R-project.org/

2.

Montoye HJ, Washburn R, Servais S, Ertl A, Webster JG, Nagle FJ. Estimation of energy expenditure by a portable accelerometer. Medicine and science in sports and exercise. 1983;15(5):403-407.

3.

Migueles JH, Aadland E, Andersen LB, et al. GRANADA consensus on analytical approaches to assess associations with accelerometer-determined physical behaviours (physical activity, sedentary behaviour and sleep) in epidemiological studies. British journal of sports medicine. 2022;56(7):376-384.

4.

Neishabouri A, Nguyen J, Samuelsson J, et al. Quantification of acceleration as activity counts in ActiGraph wearable. Scientific Reports. 2022;12(1):11958.

5.

Choi L, Liu Z, Matthews CE, Buchowski MS. Validation of accelerometer wear and nonwear time classification algorithm. Medicine and science in sports and exercise. 2011;43(2):357.

6.

Montoye AH, Clevenger KA, Pfeiffer KA, et al. Development of cut-points for determining activity intensity from a wrist-worn ActiGraph accelerometer in free-living adults. Journal of Sports Sciences. 2020;38(22):2569-2578.

7.

Van den Boogaart KG, Tolosana-Delgado R. “Compositions”: A unified r package to analyze compositional data. Computers & Geosciences. 2008;34(4):320-338.

8.

Templ M, Hron K, Filzmoser P. robCompositions: An r-package for robust statistical analysis of compositional data. Compositional data analysis: Theory and applications. Published online 2011:341-355.